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We present a method to impose linear shear flow in discrete-velocity kinetic models of hydro- 
dynamics through the use of sliding periodic boundary conditions. Our method is derived by an 
explicit coarse-graining of the Lees-Edwards boundary conditions for Couette flow in molecular dy- 
namics, followed by a projection of the resulting equations onto the subspace spanned by the discrete 
velocities of the lattice Boltzmann method. The boundary conditions are obtained without resort 
to perturbative expansions or modifications of the discrete velocity equilibria, allowing our method 
to be applied to a wide class of lattice Boltzmann models. Our numerical results for the sheared 
hydrodynamics of a one-component isothermal fluid show excellent agreement with analytical re- 
sults, while for a two-component fluid the results show a clear improvement over previous methods 
for introducing Lees-Edwards boundary conditions into lattice Boltzmann. Using our method, we 
obtain a dynamical steady state in a sheared spinodally decomposing two-dimensional fluid, under 
conditions where previous methods give spurious finite-size artifacts. 

PACS numbers: 

I. INTRODUCTION 

Planar Couette flow, a volume-preserving flow with a constant velocity gradient, is frequently used to probe the 
response of soft matter to shear [jjj. Soft materials display a wide range of interesting behaviour in shear flow: shear 
suppresses the critical temperature in a critical fluid 0; destroys long-range order in certainphases Q while inducing 
it in others and can lead to complex rheological behaviour in liquid-crystalline phases 

Couette flow is realised in the laboratory by confining fluid between two concentric cylinders which are then set 
into relative motion. Forces are transmitted from the cylinders to the fluid, setting up flows adjacent to the plates, 
which finally develop to give a linear velocity profile between them. The profile deviates from linearity in a narrow 
region near the plates called the Knudsen layer, where the fluid velocity changes rapidly. In the majority of cases, 
the width of the Knudsen layer is negligible compared to the separation between the plates and the velocity profile 
may be assumed to be linear throughout. The resulting rheological response is that of the bulk and has negligible 
contributions from boundary effects. When the radius of the cylinders is large compared to their separation, the flow 
approximates that between two parallel plates, a geometry known as planar Couette flow. 

A direct replication of such a setup in a computer simulation provides less satisfactory results: the width of the 
boundary layer or of the inhomogeneities introduced by the walls is often comparable to the size of the central layer 
where the flow profile is linear. A simulation where bulk behaviour is desired is then largely contaminated with 
boundary effects. Reducing such boundary effects by increasing the overall size of the simulation is inefficient, and 
when computational power is limited, often impossible. 

In simulations without shear flow, boundary effects are eliminated by the use of the Born-von Karman periodic 
boundary conditions (PBC). In 1972 Lees and Edwards proposed a remarkable algorithm which consists of a simple 
modification of PBC that enables the simulation of Couette flow without the introduction of physical boundaries. 
The flow has no boundary layer and no inhomogeneities are introduced. The algorithm is extremely useful when 
steady-state bulk behaviour needs to be probed in an accurate and efficient way 7]. For reasons explained below, we 
shall refer to the algorithm of Lees and Edwards as sliding periodic boundary conditions (SPBC). 

The SPBC was originally formulated for use in molecular dynamics. However, Couette flow can be described in 
at least two other frameworks, those of kinetic theory and hydrodynamics. Molecular dynamics, kinetic theory and 
hydrodynamics form a hierarchy of descriptions, in which one can proceed from the most detailed level (molecular 
dynamics) to the least detailed (hydrodynamics) by a "coarse-graining" procedure. The Bogoliubov Q method of 
passing from Newton's equation of motion to the Boltzmann equation and the Chapman-Enskog expansion Q, used 
in passing from kinetic theory to hydrodynamics, are examples of such "coarse-graining" procedures. It is reasonable, 
then, to expect that algorithms developed for molecular dynamics may find use in kinetic theory, and the resulting 
kinetic method may be useful for the simulation of hydrodynamic flow. In this paper we show that SPBC in molecular 
dynamics can be exactly represented in terms of a scattering kernel in kinetic theory, and that this kinetic boundary 
condition can be used in the lattice Boltzmann equation to provide an useful method for the simulation of the 
hydrodynamics of linear shear flow which reproduces bulk behaviour accurately and efficiently. 

There have been two earlier efforts at eliminating boundary effects in LBE with shear flow. Wagner and Yeomans 
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[T(j used modified equilibria to impose a change in momentum. Subsequently, Wagner and Pagonabarraga (WP)[Tl| 
proposed a heuristic transformation rule (which required no change in the LBE equilibria) on the one-particle distri- 
butions, based on a perturbative expansion of the distributions around the shear flow steady state. The truncation 
of the perturbative expansion implicit in their work leads to accumulative numerical error at large times. Such errors 
are avoided by the algorithm presented in this work. 

The plan of the paper is as follows. In the following section we provide a general summary of kinetic methods 
like the LBE and introduce a new kinetic scheme for the Cahn-Hilliard equation. We combine these two methods 
to produce a LBM for binary fluids. In Sec lIIII which contains the main results of this paper, we derive the kinetic 
analogue of the Lees-Edwards SPBC in the continuum using the formalism of kinetic boundary conditions. The 
difficulties associated with the discretisation of these boundary conditions on the lattice are discussed and solved. In 
Sec lIVI wc present detailed numerical results to show that our algorithm compares favourably with the earlier work 
of Wagner and Pagonabarraga. We conclude with a summary of our results and directions for future work. 



II. DISCRETE VELOCITY KINETIC MODELS 

The recent resurgence of interest in discrete kinetic models |l2i IT^ has mainly been due to the discovery that 
a specially discretised Boltzmann equation can be used as a flexible Navier-Stokes solver. This lattice Boltzmann 
equation (LBE) was originally motivated by the desire to improve lattice gas models [lj, [l5| , but was later recognised 
to be a special velocity discretisation of the Boltzmann equation [lti Il7| . The LBE represents a hyperbolic set of 
partial differential equations (PDEs) with a constant, linear advection operator whereas the Navier-Stokes equations 
are parabolic with a time-dependent non-linear advection operator. An accurate numerical solution of the former is 
easier than that of the latter. The LBE numerical algorithm is explicit and involves nearest-neighbours only, making 
it heavily parallelisable. Finally, complicated boundary conditions can easily be incorporated with little additional 
effort. Due to this combination of numerical and computational advantage, the LBE is now widely considered to be 
competitive with more traditional methods of computational fluid dynamics |13| . 

The LBE has spawned a large number of variants which collectively go under the name of the lattice Boltzmann 
method (LBM). They share the basic philosophy of using a set of hyperbolic PDEs of the form 

d t fi + c 4 • V/j = -Cij(fj - f°) (1) 

(often abbreviated as kinetic equations) to solve a set of parabolic PDEs in the form of conservation laws 

d t C a + Vj a = (2) 

for the set of conserved variables C a . The connection between the two sets of equations is established by defining 
the conserved variables to be velocity moments of the "distribution" functions fi appearing in Eq^ The matrix 
dj depends on certain symmetry requirements, while the "equilibrium" distribution functions ff depend on the 
fluxes j a which, in general, are non- linear functions of the conserved variables C a . Parabolic conservation laws arise 
in a variety of physical situations where diffusion is accompanied by advection. The incompressible Navier-Stokes 
equations representing the diffusion of vorticity and its advection by the flow field, and the Cahn-Hilliard equation 
describing the dynamics of a conserved scalar, are both of this type. 

In certain instances, the kinetic equations Eq^ follow from an underlying physically meaningful kinetic description 
where the continuum analogue of the "equilibrium" distribution function /? describes true statistical mechanical 
equilibrium. This is true in the LBE, where the continuum analogue of /? is the Maxwell-Boltzmann distribution. In 
other instances, the kinetic equations do not follow from Boltzmann kinetic theory and the "equilibrium" distribution 
functions do not have the same statistical mechanical significance as in the LBE. In such circumstances, the LBM is 
best viewed as a kinetic relaxational scheme 18] with the ff dictating the fixed point of the relaxation process and 
attempts at a statistical interpretation of the kinetic equations are best avoided. To keep this distinction in mind, we 
shall refer to relaxation schemes whose equilibria do not follow from a Boltzmann equation as lattice kinetic equations 
(LKE) and reserve the terminology of LBE for kinetic equations derived from Boltzmann equilibria. This distinction 
is important when assessing the accuracy of the SPBC for the LBE and LKE methods. Below we describe in detail 
the LBM for binary fluids, which combines the LBE for solving the Navier-Stokes equations with a new LKE method 
for the Cahn-Hilliard equation which is closely related to LBM of Swift et al . 



A. LBE for isothermal Navier-Stokes 



The set of conservation laws of mass and momentum d t p + V ■ g = and dtg + V • n = contain the hydrodynamic 
description of an isothermal fluid. They correspond to the dynamics of two conserved quantities, the mass density p, 
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and the momentum density g. The momentum density is itself the mass flux g = pv, while II Q ^ is the momentum 
flux. The fluid velocity is defined by v = g/p. The momentum flux for a incompressible Newtonian fluid is given by 

n Q /3 = g a V/3 + p5 a f3 - V^aVp + VpV a ) (3) 

Together with the conservation laws, this lead to the Navier-Stokes equations of isothermal fluid flow 

p(d t v + v • Vv) = -Vp + r/V 2 v (4) 

The lattice Boltzmann equation (LBE) is derived from the discrete velocity BGK equation 

d t fi + c i -Vf i = -£ ij (f j -f°) (5) 

by a further discretisation of space and time. The discrete velocities {c^} are the nodes of a Gauss-Hermite quadrature 
[lii Il7| which ensure that the conserved moments of the distribution functions have the same values as in the 
continuum. For the isothermal LBE, these are the mass and momentum densities, 

^2 /i{l,c M } = {p, pv a } (6) 

i 

with Greek indices denoting Cartesian directions. The Boltzmann kinetic description is restricted to a dilute gas with 
an ideal equation of state p — nksT — pc 2 sl where n is the number density and c s is the isothermal sound speed. It 
is then convenient to subtract this trivial kinetic contribution to the pressure from the momentum flux tensor and 
define a deviatoric momentum flux 

S a /3 = II Q/ 3 - pc 2 s 8 a f3 — ^2 fiQia/3 (7) 

i 

where Qic,q = Ci a Cip — c 2 5 a p and I\- a p = J^i fi c iaCi/3 is the momentum flux. In the most commonly used DdQn 
models |20t l2l| containing n quadrature nodes in d dimensions, the equilibrium distribution functions are given by 



ft 



Wi 



pV a V/3Q ial 3 

C? 2ci 



(8) 



where Wi are the set of weights defining the quadrature and repeated Greek indices are summed over. This form is 
obtained by retaining the first three terms in the Hermite polynomial expansion of the Maxwell-Boltzmann distribution 
|16| . Finally, the BGK collision matrix Cij must satisfy the constraints imposed by the conservation laws and rotational 
symmetry. In the single time approximation Cij = r~ 1 5ij, the shear (77) and bulk viscosities (£) are related to the 
relaxation time r by 77 = pc? s T, £ = {2/d)pc 2 s T. For multiple relaxation time models, the above relations remain valid 
respectively, if r is replaced by the relaxation time of the shear and bulk viscous stress [22j . 

To derive a numerical scheme, the set of coupled h ype rbolic equations in Eq. [5] must be integrated. We integrate 
the discrete BGK equation along a characteristic [23ll24| for a time At to obtain, 

fi(x + CiAt,t + At)-fi(x,t)= ds Ji(x + as,t + s) (9) 

Jo 

where J,(x, Cj,t) = — ~ fj)- The integral above may be approximated to second-order accuracy using the 

trapezium rule to give 

fi{x + CiAt,t + At)- fi(x,t) = ^[J l (x + c l At,t + At) + J l {x,t)] (10) 

This represents a set of implicit equations which would normally be solved by matrix inversion. However, the trans- 
formation to a new set of distribution functions f i [24| defined by 

At 

f i (x,t) = f i (x,t)-—J i (x,t) (11) 
renders the system explicit. Solving the above for the fi in terms of the f t gives 

/,-/0 = (i + ^C)-i(/' /0) (12) 
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Now combining Ea ll0lllT and ll2l we obtain the second-order accurate explicit LBE scheme as 

/:(x + c,At, t + At) = fat) - AtC ik (l + - /°) (13) 

In the single time approximation dj = <%/r, the above reduces to the familiar LBGK form 

fa + Ci At,t + At) = fat) (14) 

Ea llll allows us to calculate the hydrodynamic moments using the new distribution functions in terms of their earlier 
definitions, Eq|S| Taking the zeroth and first moments of Ea llll and using the constraints of mass and momentum 
conservation = JiCi Q = 0, we find that mass and momentum densities are calculated as before, with fi 

replaced by f-. However, the second moment J^i JiQiap oc (S a p — pv a V(j) is not constrained to be vanish, implying 
that that deviatoric stress must be calculated using 



S a B — 



/(1 + ^) (15) 



Since the conserved moments p and pv a which determine the equilibrium distributions /P can be calculated using 
either the fa or the /{, the fa are redundant and the numerical scheme can be formulated entirely in terms of the ][. 
Notice that fa — > fa as At — » 0, implying that fa encodes discrete lattice effects. This is true, in particular, for the 
formula for the deviatoric stress fEa ll5|l which goes over to the continuum formula (Eq|7J| as Ai/r — * 0. 

Mathematicall y , t he passage to the hydrodynamic limit can be demonstrated using the Chapman-Enskog multiple 
scale expansion |l2l | or Grad's method of moment closure 0|. Physically, hydrodynamic behaviour is obtained 
when the mean free path of the gas is small compared to the length scales associated with the flow. In the Navier- 
Stokes regime, hydrodynamic relaxation takes place predominantly through the diffusion of vorticity rather than the 
propagation of sound. For a flow with typical flow velocity v, a length scale L, and a viscosity 77, the time scales 
associated with sound propagation and vorticity diffusion are L/c s and pi? jr\ respectively. This ratio, which we define 
to be the Knudsen number, can be conveniently expressed as Kn = Ma/ Re, where Ma — v/c s and Re = pvL/rj are 
the Mach and Reynolds numbers. The LBE reproduces incompressible Navier-Stokes hydrodynamics only for Kn <C 1 
and Ma < 1. 



B. LKE for Cahn-Hilliard 



The Cahn-Hilliard equation (also known as Model B in the Hohenberg-Halperin classification of critical dynamics 

S3), 



d t tp + V ■ (V>v - MVfx) = 



(16) 



corresponds to the conservation law of a single conserved variable ip with a flux ipy — AfV/i. The chemical potential fi 
is obtained from the Landau- Ginzburg- Wilson functional F[ip] controlling the equilibrium fluctuations in ip through 
fi = SF[tp]/Sip. The mobility M is assumed to be constant and independent of tp, though the velocity v is in general 
both space and time dependent. A kinetic relaxation scheme for the above is obtained by introducing distributions 
gi obeying the kinetic relaxation equation 



d t gi + Ci ■ Vgi = -Cl {g 3 - 5?) 



(17) 



with ip = gi, and j a = 9i°ia- The model is fully determined by specifying the collision matrix , the stationary 
function g®, and the nodes of the quadrature {c-i}. In what follows we used the standard DtfQn quadrature of the 
LBE, a single-time collision matrix, Cf- = Sij /t^ and the following stationary distributions 



9i 



Wi 



(lpVaV/3 + p.S af 3)Q laf 3 



2cf 



+ 5m4> 



(18) 



where 5io = 1 for the zero- velocity population but zero otherwise. The motivation for this choice is provided elsewhere 
[26J, but note that it satisfies the following constraints 0] on its first three moments: J2i 9i = Si 9i°ia = i>v a , and 
J2i 9i c iaCip = p,S a /3 + ipVaVp. A fully discrete second-order numerical scheme is obtained by the characteristics-based 



5 



integration method introduced in the previous section. As before, this leads to an implicit set of equations for the gi, 
which can be rendered explicit by introducing the new distributions g\ 



At ^ 
9i = 9t - ~2 Ji 

where jf = —Cfjgi — g®). This gives the fully discrete LKE 



&(x + aAt, t + At)= 5 ;(x, t) - (g- - gf)/( H + ^) 



(19) 



(20) 



Using the constraint of order parameter conservation jf = we find that the order parameter is given as in the 
continuum by ip — g[, but the flux has discrete lattice corrections and must be calculated using 



3c 



\- Ai 



/(I 



At 

2t, 



-) 



(21) 



As with the deviatoric stress in Navier-Stokes fEa ll5[) the flux formula reduces to its continuum definition as Ai/r^ — > 
0. The velocity and chemical potential in the stationary distribution gf are prescribed and so do not need to be 
determined through the moments of gi. 

The formulation presented above is closely related to the binary fluid LBM of Swift et al 19]. Our model differs 
from theirs in two important ways. First, our model is formulated in the continuum from which we derive a second- 
order accurate discrete model using the method of characteristics; Swift et al begin with a discrete model which is 
only first-order accurate. Second, our stationary distributions are distinct from their "equilibrium" distributions, 
though both give identical moments up to the second. We find that with our choice of g® both isotropy and numerical 
stability is greatly enhanced. A detailed comparison of the two methods is presented elsewhere |26|; here we only 
present comparative results for shear flow using SPBC. 



C. LBM for binary fluids 



The equations of binary fluid hydrodynamics result from combining the Cahn-Hilliard equation with the Navier- 
Stokes equation with an additional force density representing the exchange of momentum between the two components 
of the fluid. We present them for completeness: 

p(d t -v + v • Vv) = -Vp + ?]V 2 v + ipVp (22) 



d t ip + V • (^v - MV/i) = 



(23) 



The order parameter ip now represents the local concentration, whose equilibrium fluctuations are governed by F[ip] = 
J cPx-lA^+Bip^+KCVip) 2 ] [23. The additional ipVp, term in the momentum equation can be written as the divergence 
of a stress tensor cr^, such that V^crTg = f/;V a /i. Then, there are at least two ways of incorporating this information 
into the LBE. The first, followed by Swift et al simply modifies the definition of the equilibrium deviatoric stress 
to include this term. Thus the equilibrium distributions for the fi are now given by 



a/3 



fi 



pVaCi a 

(pv a vp + a^ p )Q K \ ■ ) 



p + 



(24) 



The above "equilibrium" distribution can no longer be derived using Gauss- Hermite quadrature of a stationary solution 
of an underlying Boltzmann equation and is best thought of as an ansatz to be justified post-facto. Further, this 
definition breaks Galilean invariance and is responsible for the some of the artifacts seen in the model |2^,|2^. An 
alternative and more attractive method _29] is to treat the tp^P term as an external force density in the forced 
Boltzmann equation. This requires no change to the LBE equilibria and is Galilean invariant by construction. In 
this paper we follow the method of Swift et al so that a like-for-like comparison can be made with earlier shear flow 
algorithms [ll|. No change is required in the LKE for Cahn-Hilliard equation, other than identifying the velocity v 
with that in the Navier-Stokes equation. This completes our presentation of the LBM for binary fluids. In the next 
section we show how SPBC is implemented in the model which can then be used to study the hydrodynamics of a 
sheared binary fluid. 
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III. BOUNDARY CONDITIONS FOR PLANAR COUETTE FLOW 



A. Molecular dynamics 



Let us briefly recall the SPBC proposed by Lees and Edwards |6j, assuming the simulation unit cell to be a 
parallelepiped with faces at x — ±L, y = ±L and z = ±L. Periodic boundary conditions provide a prescription 
to deal with coordinate values that flow out of this range as the integration of the equations of motion proceeds. 
PBC simply identifies all coordinates modulo the size of the unit cell 2L: any coordinate which lies outside the range 
[— L, L] is mapped back into it using modular division. Additionally, the PBC is applied to each coordinate separately 
and does not affect the velocities of the particles. Thus, with PBC applied in the y direction, a particle leaving the 
top face at y = L is simply reinserted at the bottom face at y = —L, with no change to the y and z coordinates or 
the velocity. Similar rules apply for PBC in the x and z directions. 

The SPBC proposed by Lees and Edwards follows from a simple modification of PBC. It prescribes that a particle 
leaving a face of the box in the velocity-gradient direction is re-entered at the opposite face but at a new position 
displaced along the flow direction, and with a new velocity incremented only in the flow direction. The ordinary 
PBC is applied for particles leaving the box in the remaining directions. For concreteness, assume that the flow 
and velocity-gradient directions are x and y respectively. The flow velocities at ±j/ = L are ±t/, giving a velocity 
gradient 7 = U/L. The SPBC then requires that a particle leaving the top face at y = L with coordinates (x, L, z) 
and velocities (c x , c y ,c z ) is reinserted at the bottom face y = —L with a new x-coordinate x' and new x-component 
of velocity c' x given by 



so that its final coordinates and velocities are (a;', —L, z) and (c' x , c y , c z ) respectively. Particles leaving the box through 
the bottom face at y = —L re-enter at the top face at y — L using the same rule as in Eq|23but with the sign of U 
reversed. The standard PBC is applied to particles leaving the box through the faces at x = ±L and z = ±L. 

The motivation for these rules becomes transparent when we consider the deformation of a cube with edges of 
length L in uniform shear flow with 7 = U/L. The cube undergoes an affinc deformation where the top and bottom 
faces move in opposite directions with relative velocity 2U . Points at the top and bottom faces which had identical 
x-coordinates at t = are displaced by 2Ut at time t. 



The Lees-Edwards sliding periodic boundary conditions described above operate on the coordinates and velocities 
of the particles crossing the boundary dV of the simulation volume V. The kinetic description of a molecular system, 
on the other hand, is formulated not in terms of particle coordinates and velocities, but in terms of the one-particle 
distribution function /(x, c,i), representing the number density of particles at the phase point (x, c) at time t. 
Accordingly, to make use of SPBC in kinetic theory, we must translate the SPBC condition Eq|23into an equivalent 
condition on the one-particle distribution function. This can be done using the scattering kernel approach for reflective 
boundary conditions in kinetic theory |3C| . as we show below. 

Consider, as before, fluid flow in a domain V whose boundary is dV . Let s 6 dV be a point on the boundary 
and n(s) be the inward pointing normal at that point. Molecules at the phase point (s,c) leave V and impinge 
on the boundary if n(s) • c < 0; similarly, molecules enter V and emerge from the boundary if n(s) • c > 0. Let 
B(s,c — > s',c';t) be the probability that at time t, a molecule which strikes the boundary at point s with velocity c 
emerges at point s' with velocity c'. For a non- absorbing wall, a molecule is always re-emitted, which provides the 
normalisation condition for P 



The total number of molecules emerging from the boundary (per unit time and area) at the point s' with velocity 
c' must be equal to the number of molecules of all velocities incident on all parts of the boundary that change their 
velocity to c' and emerge at the point s'. Stated mathematically |30j . 



x' = x - 2Ut, 4 



(25) 



B. Kinetic theory 




(26) 




dsdcB(s, c — » s', c'; t)Q(—c ■ n)|c ■ n|/(s, c, t). 



(27) 



we obtain the boundary condition on the populations emerging from the boundary in terms of the incident populations 
and the scattering kernel B which characterises the particle- wall interaction. The Heaviside step function selects the 
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impinging and emerging populations from the total populations at the boundary. The majority of physical boundary 
conditions are local in space and independent of time, yielding a scattering kernel B that is dependent non-trivially 
on the velocities only. In the simplest case of specular reflection, the particle velocity normal to the wall is reversed 
at the point of contact. In the opposite case of diffuse scattering, the particle is thermaliscd on contact and leaves 
the wall with a velocity sampled from a Maxwellian distribution determined by the wall temperature. The scattering 
kernel for specular scattering is B(s, c — » s', c') = 5(s — s')6(c + 2n(n • c) — c') while that for diffuse scattering 
is B(s,c — > s',c') = 6 (a — s')|n • c| exp(— c 2 /fcsT)/27r(fcBT) 2 More realistic wall-particle interactions are obviously 
possible and have been explored in detail |30| and implemented in the LBE |3lj . 

Viewed in this light, Lees-Edwards boundary conditions can be thought of as a non-local scattering process, in 
which particles incident on the upper boundary cross over to be re-emitted from the lower boundary, with the initial 
and final coordinates and momenta being given by Eg 1251 It is clear, then, that the transformation rules for the 
particle coordinates and velocities at the boundary can be used to construct the scattering kernel B. Enforcing the 
SPBC rules fEci l25fl at the boundaries through delta functions, we find that 

B S pbc(x, c -» x', c'; t) = S(x ± 2Ut - x')5(y + y')S(z - z')S(c x ±2U - c' x )6(c y - c' y )5(c z - c' z ) (28) 

Inserting this into the Eg 1271 and completing the integration, we obtain the equivalent of Lees-Edwards boundary 
conditions for the distribution function: 

f(x, y, z, c Xl c y , c z ,t) = f(x ± 2Ut, -y, z, c x ± 2U, c y , c z ,t) (29) 

Below we shall implement this SPBC on the distribution functions of the LBE and LKE, paying special attention to 
the difficulties in going from a space where x, c are continuous to a discrete velocity space and a discrete space-time 
lattice. 



C. Discrete implementation 

The SPBC, EqEHI poses no problem so far as the temporal discretisation is concerned, since it relates values of the 
distribution function at equal times. However, a direct implementation of Eo 1291 on the lattice is complicated by the 
fact that the nodes of the quadrature (ie. the set of velocities Cj) are discrete, as is the position x. The discreteness of 
the velocity space must be dealt with even if the model is formulated in the continuum, as is the discrete velocity BGK 
equation Eq|S] For such a model to be computationally useful, it must ultimately live on the lattice, requiring us to 
address the discretisation of space as well. Thus there are two distinct, and conceptually separate problems: we deal 
with the velocity discretisation first. Let us first recall that any (well-behaved) function can in general be expanded 
in terms of orthogonal polynomials. A particularly useful expansion for the / is in terms of Hermite polynomials |32j |. 
where the coefficients of the expansions are Hermite moments of the distribution function. This expansion is the basis 
of the derivation of the discrete velocity BGK equation from the Boltzmann-BGK equation For an isothermal 
LBE, at least three terms must be retained in this expansion. The expansion coefficients are evaluated by Gaussian 
quadrature 01 > the nodes of the quadrature defining the discrete velocities Cj . The (x, t) can then be represented 
in terms of the velocity moments of /(x, c, t) as 



ft = m 



pV a Ci a S a pQi a p 



(30) 



where Gi represent higher moments of the distribution function, variously known as kinetic modes, non-hydrodynamic 
modes, or "ghost" modes [33|, |34|, |35| . We circumvent the problem of the discreteness of velocity space as follows: 
instead of implementing the SPBC Eq[2Hl directly in terms of the distributions /(x, c, t), we obtain equivalent condi- 
tions on their moments, which are then "projected" on the discrete velocity space by using Ea l3UI Taking the first 
three moments of Ecil29lfollowed by a change of variables in the integrals we obtain the SPBC conditions on the mass, 
momentum and stress densities as: 

P = p' (31) 
pv a = (pv a )'±p'U a (32) 
S af3 = S' a p±(pv a )'U l3 ±(pv f3 yU a + p'U a U( } (33) 

where the p — p(x, y, z, t), p' = p{x ± Ut, ^fy, z, t) etc, and U a — 5 xa U. We now demand that the first three moments 
of the Eq|2Hbe satisfied at the discrete level. This gives us the discrete analogue of the SPBC for the LBE. 

fi(p,pv a ,S a p,G) = .n(p' : (pv a y,S' af3 ,G) (34) 



8 



The transformation of the model dependent kinetic modes G; have been ignored for clarity of p resentation, but can 
easily be included in the above framework. (We list all the modes for the DdQ15 model in |29|: for DdQ9 see [34|). 
Below we discuss why neglecting the transformation of the kinetic modes may be justified. 

The SPBC for the LKE for the Cahn-Hilliard equation can be implemented by repeating the steps above. First we 
note that the gt can be written as 



<], 



at 

2cf 



Gj 



Sioip 



(35) 



where /Lt Q| g = '^2 i giQiap, and Gf represent the kinetic modes. Since the transformation rule for the moments of 
the distribution function is determined by completely by the number of velocities that enter into its definition, 
the transformation rules for the LKE moments follow identically from the counterparts in the LBE by making the 
replacements p — > ip, (pv a ) ja, and S a p — > fi a p- We thus have: 



V> = ip' 

30i = j'a ± i>'U a 
Ma/3 = Ma/3 ± j'oPp ± j'/3 Ua + *P'U a Up 

Leaving the kinetic modes untransformed, we obtain the discrete SPBC for the LKE as 

giilpJa^ap,^) = g'0 ,3'^^^) 



(36) 
(37) 
(38) 



(39) 



A linear interpolation is applied to deal with the case when x — Ut is not a member of the lattice. The distribution 
function at the point x — Ut is partitioned between the nodes bracketing the point x — Ut in proportion to their 
distance from the point. With the transformation rules and the linear interpolation we now have a complete discrete 
SPBC algorithm. Notice that we arrive at the algorithm from the continuum SPBC using two approximations: 
first, we do not transform the kinetic modes and second, we apply a linear interpolation to deal with the spatial 
discreteness. The first of the approximations is easily avoided, once expressions for kinetic modes are known in terms 
of the microscopic velocities. In the hydrodynamic limit, the kinetic modes are only weakly coupled to the conserved 
variables |3(J and their transformation may be neglected to a first approximation. In the interests of simplicity, we 
have only presented the model- independent transformation rules for the hydrodynamic modes. The second of the 
approximations is unavoidable on a discrete lattice, though more accurate interpolation schemes might be envisaged. 
As we show in the following section, the present combination already gives results of acceptable numerical accuracy. 



D. Large shear rates 

The LBE constraint of small Mach numbers (U « c s ) places an upper limit on the shear rate 7 = U/L which 
becomes progressively smaller with increasing system size L. To overcome this limitation, Wagner and Pagonabarraga 
|10| introduced internal Lees-Edwards boundaries, where the Lees-Edwards boundary conditions are applied without 
the application of PBC. In the particle picture, this corresponds to particles crossing the boundary accompanied by 
a change in velocity ±f and a change in the x— coordinate ±f t, but no change in the y and z coordinates. In terms 
of distribution functions, we have 

f{x,y,z,c x ,c y ,c z ,t) = f{x±2Ut,y,z,c x ±2U,c y ,c z ,t) (40) 

whose discrete implementation follows as in the previous section, except that PBC is no longer required in the y— 
direction. With the introduction of Nle such planes, the maximum attainable shear rate is Nle times that with the 
SPBC. In the following section we demonstrate the use of LE planes in the study of spinodal decomposition under 
shear. 



IV. TEST APPLICATIONS 



We now provide results of numerical simulations using the SPBC derived in the previous section, treating in turn 
the hydrodynamics of sheared one-component and two-component fluids. In what follows, we define lattice units of 
time, length and mass by setting At = 1, Ax = 1, po = 1, where po is the mean mass density at a lattice node. 
In our simulations, we use a collision operator which relaxes the stress modes at a rate r _1 and projects out all the 
non- hydrodynamic modes |37l ]. p2l |. For a simple fluid the only adjustable parameter, then, is the relaxation time r, 
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now related to the shear viscosity by rj = Pqc 2 t. For a binary fluid, the adjustable parameters are the mobility M, the 
interfacial thickness £ and the interfacial surface tension a. In the LKE the mobility is related to by M = c 2 t^, 
while the parameters contained in the LGW functional F[ip} determine £ and a [221]. Our choice for all the binary 
fluid parameters is that of Kcndon et al 27]. In all cases, the simulation is performed in a rectangular domain of 
dimension L x x L y . Shear is applied so that flow takes place along the x— axis and the flow gradient is along the 
y— axis. 



A. Single fluid 



For the single fluid we validate our method by comparing against two cases where analytical solutions are known. 
The first of these is impulsively started Couette flow, where a quiescent fluid is suddenly subjected to shearing motion. 
The relaxation of the velocity to its asymptotic linear shear profile is determined by the diffusion of momentum from 
the boundary into the bulk: the solution for flow between parallel plates moving at relative velocity U with separation 

l is m 

2u ( 2 2 v ^ \ ■ n7T y 



=i 



v x (x, y, z, t) = U(l - y/L) - — ^ exp I -nV — \ sin (41) 



where the kinematic viscosity v = rj/po is identical to rj in our choice of units. The SPBC algorithm can be used to 
study transient phenomena if the time scale of the transient is large compared to the time it takes for vorticity to 
diffuse across the transverse dimension of the system. To compare this with the results of an LBE simulation we must 
ensure that that the Knudsen number Kn = v/Lc s is small, so that the the vorticity diffusion time scale L 2 jv is much 
larger than the sound propagation time scale L/c s , and incompress ible hydrod yna mic b ehaviour is obtained. Three 
representative results from a wide range of viscosities shown in Figs 1 (a) l(b)| and 1 (c) where the numerical results 



are compared with the theoretical expression in Eg 1411 Table [I] shows the relative error e v = (v n — v t )/v t , where v n is 
the numerical results, and v t the theoretical expression, half a lattice site away from the LE plane. As can be seen, 
the error increases substantially with decreasing viscosity and decreases with time as the velocity profile relaxes to its 
steady state. Overall, errors remain very low in most situations apart from the very early stages of the most inertial 
runs. For this test, we find very little difference between our algorithm and that of the Wagner and Pagonabarraga 
for the range of parameters we have have explored; accordingly we do not display data for their algorithm. 

Our second test explores the effect of a Gallilean transformation along the direction of the velocity gradient, or 
equivalently, that of an initial condition in which the fluid is moving uniformly with a velocity U 1 - in the direction 
of the velocity gradient. In this case, we expect that the momentum of the system in the flow direction will increase 
linearly with time. This is most clearly understood in terms of SPBC for molecular dynamics: the boundary condition 
does work on each particle as it crosses the boundary, since it changes the velocity of particle in the flow direction 
by ±U. A useful description is in terms of Lees-Edwards "demons" sitting at the boundary and giving a transverse 
impulse to every particle crossing the boundary. Since the total number of particles in the simulation box is constant, 
the number of particles leaving from the boundary is exactly balanced by the number of particles entering through the 
boundary. This implies that summed over the incoming and outgoing particles, the demons add equal but opposite 
amounts of momentum to the system, resulting in no change of momentum in the flow direction. (However, the kinetic 
energy, being proportional to U 2 does not cancel; the demons add energy to the system which is dissipated as heat.) 
However, with a net transverse flow, there is a net flux of particle crossing the boundary in a given direction and the 
momentum added by the demons no longer cancels out, resulting in an increase (or decrease) of momentum along the 
flow direction. For a constant transverse flux, the same amount of momentum is added at every time step, resulting 
in a linear increase of the transverse momentum. The same argument applies in kinetic theory. In the LBE, the 
energy density is not conserved, so the effects of heating are not noticeable. However, the argument for the increase 
in momentum still holds. Formalising our previous argument we see that the total momentum should increase as 

f'i (x, t)cte = P * = PoU^L x Ut (42) 

x,i x 

which is simply the product of the number of particle crossing the boundary poU^Lx, the momentum added per 
particle U, and the duration of this process t. This should lead to a steady drift with time of the entire velocity 
profile. 



Our results for this test are shown in Fig 2(a) The system was initialised with a linear velocity profile with a 
transverse velocity U± and all three components of the total momentum was measured. As expected, there was no 
change in the or z— components, while the x— component increased linearly with time with a slope predicted by 
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Eg 1421 Further, the entire velocity profile shifted by at a constant rate as shown in Fig 2(b) This validates the 
Gallilean invariance of the method. Similar simulations performed with the algorithm of Wagner and Pagonabarraga 
show no increase (the a;— momentum remains constant) which clearly demonstrates the lack of Gallilean invariance. 
This has been noted in previous work where a sheared binary fluid simulation with a net transverse momentum 
produced unphysical results |39| . We have verified that mass was properly conserved in both our method and that of 
Wagner and Pagonabarraga. 



B. Binary fluids 



For the binary fluid, we have performed tests for steady-state deformation of a droplet in shear flow and spinodal 
decomposition under shear. In either case, to the best of our knowledge, no analytical results are available in the 
range of parameters accessible in a LBM simulation. Therefore, we must satisfy ourselves with qualitative tests based 
on a visual inspection of the data, focusing on possible artifacts near lattice planes where SPBC is applied. In the 
binary fluid case a like-for-like comparison with WP requires care since our binary fluid LBM differs significantly from 
theirs, even before the different handling of SPBC is introduced. Our preliminary investigations indicate that the g® 
used by WP introduce the largest source of error over and above those that may be introduced by their boundary 
conditions. Accordingly, we present results for spinodal decomposition and droplet deformation for the WP algorithm 
using the LBM outlined in their original work (henceforth WPO), the WP algorithm using the stationary distributions 
in EqHH| (henceforth WP1) and our LBM algorithm with SPBC as presented in SecE] (henceforth SPBC). 

For our first test we looked at the steady-state deformation of a sheared droplet comparing SPBC with WPO. For 
a droplet of radius R and surface tension er in a shear flow with shear rate 7 the important dimcnsionlcss groups are 
the Reynolds number at the scale of the droplet Re = p'jR 2 ji) and the capillary number Ca — rjjR/cr. The capillary 
number is the ratio of the viscous stress to the surface tension. The system was sheared using one SPBC plane with 
U = 0.032 for system of size 128 x 128, resulting in a shear rate of 7 = 2.50 x 10 -4 . In the initial state a droplet of 
radius R = 32 was centred respect to the SPBC boundary. The fluid viscosity was chosen such that Ca = 0.038 and 
Reynolds number Re — 1.28 with their ratio being Re/Ca — 33.6. The droplet was then visualised in the steady-state 
(t = 100000). Fig. clearly shows the existence of artifacts where the LE plane crosses the droplet. In particular, 
one notices the presence of a cusp in the visualisation in both WPO and SPBC schemes, but this artifact significantly 
more marked with WPO than with our new scheme. In the isoline plot of Fig. 01 spurious behaviour is clearly observed 
at the boundary in the WPO scheme, but is almost unobservable SPBC. 

Our next test involves spinodal decomposition in shear flow. We performed two sets of simulations on a 320x512 
lattice with 8 and 32 LE planes. In either case, the shear rate was fixed at 7 = 0.001, requiring U = 0.04 for 8 planes 
and U = 0.01 for 32 planes. Initial conditions were set for a 50/50 mixture, with a random composition at each 
site while fi distributions were initialised to their steady state values. The remaining parameters were 77 — 0.025, 
A = -B = 0.00625, k = 0.004, M = 4.0 and a = 0.0042. Physically, the results should depend only on 7 and 
not on the number of LE planes. Pseudo-colour plots of the domain morphology at late times are shown in Fig[5] 
(SPBC), Fig© (WPO), and FigH(WPl). The morphology in Fig© and FigEJdepends on the number of LE planes. 
In particular, increasing the number of LE planes increases the number of domains that have "wrapped" around the 
flow direction. We believe this arises due to the error introduced at the LE planes by previous algorithms, as noted 
in |39j. The nature of this error is to introduce a spurious alignment of the interface with the LE planes, leading to a 
pinning of the domains and so promoting the wraparound along the flow direction. This pinning effect is clearly seen 
in Fig|S|for the WPO and WP1 methods, but is not discernable using SPBC. 

A more quantitative test follows from length scale measures L\ and Li extracted from the gradients of the order 
parameter field 0,E3- Physically, these denote the principal axes of spherical droplet deformed by the flow, reflecting 
the effect of elongation (which stretches the droplet into an ellipse whose major axis is aligned with the flow) and 
rotation (which rotates the principal axes of the ellipse). By convention, L\ is the minor axis and Li the major axis. 
Applied to a spinodal pattern, L\ measures the deformation approximately normal to the flow, while Li measures the 
deformation approximately in the direction of flow. Simulations in which the domains wrap around the flow direction 
( L<z > L x ) are contaminated by finite-size effects; no quantitative conclusions can be drawn from such data. In FigEI 
FigllOland Fig^Jwe compare the length scales L\ and Li for SPBC, WPO, and WP1 respectively. In the latter two 
cases, Li gradually increases beyond L Xl indicating wrap around, and no steady state is reached. Further, the rate 
at which the wraparound proceeds (the slope of the Li(t) curve) increases with the number of LE planes. In the case 
of SPBC we obtain a steady-state with Li < L x , with no dependence on the number of LE planes. This is the first 
observation, within an LBM simulation using periodic boundary conditions, of a dynamical steady state in which the 
domain scales L\ and Li achieve saturation in a manner that is not affected by finite size artifacts and/or lock-in 
between fluid interfaces and Lees-Edwards planes, as visible in Fig[S] More detailed results for this problem will be 
presented elsewhere [il) . 
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c. 



Limitations 



Combining these results, we feel confident that the errors introduced solely due to SPBC and internal LE boundaries 
are now under good control, and no larger than other sources of error [27j in the LBM for binary fluids. This could 
not be said for the previous state of the art as represented by WPO. However, shear flow may introduce errors beyond 
those mentioned in [23. The LBM model described here is not Gallilean invariant. This deficiency is likely to show 
up more acutely in the presence of shear flows, in particular with multiple LE planes. The small artifact at the SPBC 
planes noted in our droplet test may be due to the lack of Galilean invariance. Our SPBC rules do not break Galilean 
invariance, however. The interpolation of the distribution functions, while respecting the conservation laws, may 
introduce numerical errors. We have chosen the simplest linear interpolation scheme, and no doubt more accurate 
interpolations schemes can be used to improve accuracy. Finally, in shear flow there is always a mean advection 
velocity, implying that errors introduced by the advection scheme will be accumulative. The advection scheme of the 
LKE is first order and has a diffusive error which depends on the Courant number C = UAt/ Ax. While a systematic 
error analysis remains to be done, we have noted that acceptable accuracy is obtained only with C < 0.05. For larger 
values, the simulation is contaminated by the diffusive error of the advection scheme. 



To summarise, we have presented a new scheme for the shear flow in lattice Boltzmann and lattice kinetic models, 
which follows from a systematic coarse-graining of the Lees-Edwards boundary conditions of molecular dynamics. The 
boundary conditions are then implemented on the discrete kinetic space. The problem of discrete velocities is solved by 
applying the transformation to the moments of the distribution functions, while the problem of discrete space is solved 
by a numerical interpolation. We have validated our method for one-component and two-component fluids and find 
excellent agreement with analytical results when available, and qualitatively correct behaviour in more complicated 
cases. The present algorithm allows for the investigation of the effects of shear on a wide variety of complex fluid 
flows. The present method may be combined with our fluctuating lattice Boltzmann equation (FLBE) |36| to study 
the interplay of thermal fluctuations and shear in complex fluid hydrodynamics. We are presently investigating the 
effects of shear on spinodal decomposition in binary fluids |4l| . 



We thank Ignacio Pagonabarraga, and Paul Stansell for useful discussions and Mike Cates for helpful suggestions 
and a critical reading of the manuscript. R. A. and K. S. are supported by EPSRC GR/R67699 (Reality Grid). 
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TABLE I: Relative error in the velocity e„ (defined in the text) as a function of dimensionless time t v = r]tjl? y , half a lattice 
site from the SPBC plane for different values of the viscosity n. 
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(c) ri = 0.0005 



FIG. 1: Relaxation to steady-state in impulsively started planar Couette flow using the SPBC algorithm. The horizontal 
axis represents the y axis for the lattice, the vertical axis corresponds to the flow velocity v x (y) — J2 X v x (x, y)/L x . Curves 
within each set correspond to different dimensionless times t v — r/t/Ly. The symbols represent numerical data, dotted lines 
the theoretical expression. 
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FIG. 2: Gallilean invariance of the SPBC algorithm. (a)Variation of the total fluid momentum in the flow direction as a 
function of time with flow in the velocity-gradient direction. Data is shown for U = {0.01; 0.02; 0.04} and U = {0.015; 0.03}). 
The vertical axis corresponds to ^2 x {pv x )/ poU ± L y . Each line corresponds to the predicted values obtained using eqn l41l where 
each line corresponds to runs with identical U. Note that symbols are overlap for any given U and time, (b) Velocity profile 
v y — f(x) with U = 0.04 and U ± = 0.004. The horizontal axis represents the flow direction, the vertical axis represents the 
average velocity component v x at a location y. Simulation data has been plotted for At v 




(a) 



(b) 



FIG. 3: Surface plot z = xj)(x,y) of the steady-state order concentration ip under shear at t = 100,000 for (a) SPBC and (b) 
WP0. Notice the sharp discontinuity in the order parameter field at the SPBC plane in (b). 




FIG. 4: Steady state droplet deformation at t = 100,000 for (a) SPBC and (b) WPO. Isolines start for $ = -0.95 with steps 
—0.01 inside the droplet, and for $ = 0.95 with steps 0.01 outside. Artifacts at the SPBC plane indicated by the dotted line 
are clearly visible in the (b). 
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FIG. 7: Bicontinuous domains at step 90,000 for WP1 with (a) 8 planes and (b) 32 planes. Notice that domains are wrapped. 
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FIG. 9: Length scales (a) LI and (b) L2 for SPBC, showing no discernable dependence on the number of planes. 



70 
60 
50 
40 
30 
20 
10 










8 LE plane 
32 LE plane 


s 

s 































































1000 



800 



600 



400 



200 





I 

3 


8 LE plafiet 
2 LE plahe 






























i 







20000 40000 



60000 80000 



100000 



20000 40000 60000 80000 100000 



(a) 



(b) 



FIG. 10: Length scales (a) L\ (b) L2 for WP0 showing that L2 is weakly dependent on the number of planes. 
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FIG. 11: Length scales (a) L\ and (b)Z/2) for WP1 showing that L2 depends strongly on the number of planes. 



